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Motivated by the problem of the evolution of bulk gravitational waves in Randall-Sundrum cos- 
mology, we develop a characteristic numerical scheme to solve 1+1 dimensional wave equations in 
the presence of a moving timelike boundary. The scheme exhibits quadratic convergence, is capable 
of handling arbitrary brane trajectories, and is easily extendible to non-AdS bulk geometries. We 
use our method to contrast two different prescriptions for the bulk fluctuation initial conditions 
found in the literature; namely, those of Hiramatsu et al. and Ichiki and Nakamura. We find that if 
the initial data surface is set far enough in the past, the late time waveform on the brane is insensi- 
tive to the choice between the two possibilities; and we present numeric and analytic evidence that 
this phenomenon generalizes to more generic initial data. Observationally, the main consequence 
of this work is to re-affirm previous claims that the stochastic gravitational wave spectrum is pre- 
dominantly flat Ogw oc /°, in contradiction with naive predictions from the effective 4-dimensional 
theory. Furthermore, this flat spectrum result is predicted to be robust against uncertainties in 
(or modifications of) the bulk initial data, provided that the energy scale of brane inflation is high 
enough. 
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I. INTRODUCTION 

It is well known that the Randall-Sundrum (RS) 
braneworld model 0,0 is in excellent agreement with 
general relativity at low energies. This is the principal 
appeal of the model; it is one of the only examples of 
a scenario involving a large extra dimension that entails 
no serious conflicts with general relativity. However, this 
means that one needs to consider high energy or strong 
gravity scenarios to properly test the model. One possi- 
bility is to examine the high energy epoch of braneworld 
cosmology, where exact solutions of the 5-dimensional 
field equations are known. Well understood braneworld 
phenomena 0, review] include a modified cosmic expan- 
sion and early times and 'dark radiation' effects, whereby 
the Weyl curvature of the bulk projected on the brane 
acts as an additional geometric source in the Friedmann 
equation. 

But if one wants to move beyond the exact descrip- 
tion of the background geometry in these cosmological 
models, there are significant technical difficulties. A cos- 
mological brane is essentially a moving boundary in a 
static 5-dimensional background — anti-de Sitter space 
in the RS model — so perturbations are described by 
bulk wave equations with boundary conditions enforced 
on a non-trivial timelike surface. While it is possible to 
make some analytic progress when the brane is moving 
'slowly' [lUliH, ^ e morc interesting case of a fast- 
moving, high-energy brane remains impervious to such 
treatment. 

The purpose of this paper is to present a new numeric 
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algorithm to solve wave equations in the presence of a 
moving boundary. For the sake of simplicity, we restrict 
ourselves to a class of wave equations and boundary con- 
ditions that correspond to tensor, or gravitational wave 
(GW), perturbations. This is not the first attempt to 
deal with these equations numerically: previous efforts 
include pseudo-spectral 0, El and direct evolution 
[Til Il2l IT3I Hil ITa ] methods using various null and non- 
null coordinate systems in which the brane is stationary. 

So why introduce another method? Our primary mo- 
tivation is to develop an algorithm that offers several im- 
provements to the preexisting efforts. Our technique is 
based on the highly successful characteristic integration 
methods from black hole perturbation theory . These 
offer several advantages, not the least of which is an ele- 
gance of implementation that takes the causal structure 
of the spacetime explicitly into account. We also work 
in Poincare coordinates, which greatly simplify the bulk 
wave equation and avoid the coordinate singularities that 
plague Gaussian normal patches. This makes our algo- 
rithm both transparent and unique: while other groups 
carry out their analysis in Poincare coordinates, they al- 
ways transform the brane to a static location for actual 
calculations. Our procedure is developed from first prin- 
ciples, and we pay careful attention to discretization er- 
rors. Hence we have a good theoretical understanding of 
the convergence properties of our method, which can then 
be tested in actual calculations. The fact that the code 
behaves as expected — with explicit quadratic conver- 
gence — imparts a certain level of confidence in conclu- 
sions drawn from our numerical results. Finally, our tech- 
niques should be easily adaptable to other braneworld sit- 
uations; i.e., more complicated bulk geometries, sophisti- 
cated specification of initial conditions, multiple branes, 
etc. 
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Our secondary motivation stems from the fact that the 
numerical results found in the literature do not seem 
to agree with one another. In this paper, we limit the 
discussion to tensor type perturbations, so the princi- 
pal observational consequence of our work is the present- 
day spectral density Ogw of relic GWs generated dur- 
ing inflation on the brane [H E d M, El El El 
This 'stochastic GW background' is potentially observ- 
able by next-generation detectors such as the Big Bang 
Observatory. (Maggiore |24| offers comprehensive review 
of the stochastic GW background from a 4-dimensional 
perspective.) Hiramatsu et al. (henceforth HKT) 
have predicted that f^cw oc / for frequencies above 
some threshold f c , whose value is determined by the 
curvature scale of the bulk. On th e ot her hand, Ichiki 
and Nakamura [Til IT^ (henceforth llNt) have predicted 
Hqw fx /~ - 46 from their simulations using a differ- 
ent initial condi tion . Recently, Kobayashi and Tanaka 
P3 (henceforth IktIi have applied a different numerical 
method to the quantum mechanical version of the prob- 
lem, where one treats the entire evolution of tensor modes 
during inflation and radiation-domination as a particle- 
production phenomenon. They also derive an approxi- 
mat ely flat spectrum for high frequencies, in agreement 
with iHKTl 

What is the right answer? It is difficult to compare 
these calculations directly because each group uses a dif- 
ferent prescription for dealing with initial conditions and 
a different numerical method. It is sensible that when 
trying to solve a problem numerically, one should con- 
firm that several independent methods yield the same 
results under the same circumstances. Only then can 
we be confident in the predictio ns. To this end, we at- 
tempt to reproduce the results of lHKTl and|S| using our 
numerical method and their respective choices of initial 
condi tions. 1 We find that our numerics reproduce the 
HKT results within an acceptable tolerance, but we are 
unable to duplicate the GW spectrum predicted bv lINl 
Indeed, we find that as long as simulatio ns are started 
in sufficiently high energy epochs, both the HKT and IN 
initial conditions lead to the same flat spectrum. 

This is an inter esting and som ewha t unexpected result: 
Superficially, the lHKTl IinI . and lKTl formulations appear 
quite different from one another, yet they ultimately pro- 
duce the same late-time behaviour. This leads to the 
question: How robust is the prediction of a flat GW spec- 
trum to arbitrary modification of the initial conditions? 
Answering this is the third motivation for this paper, and 
is a fairly ambitious goal. This is because 'arbitrary mod- 
ification' implies the need to consider an infinite number 
of cases, which is highly impractical. Hence, we need 
to settle for qualitative conclusions drawn from numeric 



1 The Wronskian formulation favoured by iKTl is sufficiently dis- 
tinct from the other methods to defer its consideration to future 
work. 



simulations of individual cases coupled with some ap- 
proximate analytic results. 

The layout of the paper is as follows: In Sec. [H] we 
describe the problem we are going to solve in both gen- 
eral terms, and as specialized to tensor perturbations in 
RS cosmologies. The numeric algorithm used throughout 
the paper is developed in Sec. IIIII and then comprehen- 
sively tested in Sec. lIVI The issue of initial conditions for 
G Ws in b raneworld cosmology is reviewed in Sec. [3 and 
the HKT] and |S| approaches are explicitly contrasted in 
Sec. lVIl A more generic class of initial data is considered 
in Sec. IVIll using a combination of numeric and analytic 
methods. Sec. I Villi is reserved for our conclusions. Ap- 
pendix previews the jargon associated with the stochas- 
tic GW background and how to convert the results of 
numeric simulations into observational predictions. 

We use units in which h = c = 1 and the 'mostly 
positive' metric signature. 

II. STATEMENT OF THE PROBLEM 

A. Generic formulation 

In this subsection, we define the generic type of prob- 
lem that we solve in this paper. Consider the following 
wave equation: 

l-D 2 + V(z)]yj(t,z)=0. (1) 

Here, D a is a covariant derivative operator on a flat 2- 
dimensional manifold: 

ds 2 = -dt 2 + dz 2 = -du dv, (2) 

where u = t — z and v = t + z are the usual retarded 
and advanced time coordinates. Our goal is to obtain 
the value of the field throughout a finite region S7 of the 
(i, z) spacetime, which is depicted in Fig.Q] 

As can be seen in this figure, the boundary of £1 is 
composed of three distinct parts: d^l ± are null surfaces 
to the future and past, while dflb is a timelike surface 
that we will refer to as the 'brane'. The brane is defined 
parametrically via the equations 

t = t b (rj), z = z b (n). (3) 

The parameter n is selected to be affine in the flat (i, z) 
geometry: 

u ■ d = %d t + hd z , u - u = —1. (4) 

Here and henceforth, we have use an overdot to denote 
d/dr]. As long as z b ^ constant, the brane can said to 
be 'moving' with respect to the (t, z) frame of reference. 
We also define a normal vector to dQb pointing into O: 

u • d = Zbd t + ibd z , n-n = l, u • n = 0. (5) 

In order to have a well-defined hyperbolic problem, we 
must also specify initial and boundary conditions. The 
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such that 




FIG. 1: (Colour online.) The spacetime domain £1 over which 
we seek a numeric solution of the wave equation. Superim- 
posed on Q is a (particularly coarse) example of the compu- 
tational grid we use to discretize the problem 



initial conditions consist of choosing the value of the field 
ip on the past null boundary dQr . This initial profile can 
be selected arbitrarily. We take the boundary conditions 
on the brane to be 



= [(n•D)V>-a(?7)V>]&• 



(6) 



At this stage, we leave a as a finite, but otherwise ar- 
bitrary, function of 'time' on the brane. Therefore, this 
represents a wide variety of possible boundary conditions, 
except for pure Dirichlet ip(t b , z b ) — 0. 



B. Application: gravitational waves in brane 
cosmology 

In this subsection, we show how the equations gov- 
erning tensor perturbations in RS one brane cosmologies 
can be written in the general form introduced in Sec. Ill Al 
above. 

A cosmological RS brane is identified as a hypersurface 
in 5-dimension anti-deSitter space with cosmological con- 
stant -6/£ 2 : 



ds 2 = -^(-dt 2 + 5 lj dx l dx 3 + dz 2 ). 

The brane is defined by the parametric equations 
t = t b {rj), z = z b (r)), 



(7) 



(8) 



t b = Jl + z 2 



(9) 



This equation ensures that the metric on the brane, 
ds 2 = a 2 (rj)(~dri 2 + 5ijdx l dx j ), 



(10) 



is of the standard cosmological form with rj as the con- 
formal time, and the scale factor identified as 



a{rj) = l/z b (ri). 



(11) 



The brane's dynamics are described by a Friedmann 
equation derived from the Israel junction conditions: 



p(a) 



2A 



(12) 



Here, p is the density of brane matter, k 2 = SttG^ is the 
4-dimensional gravity-matter coupling, and A is the brane 
tension. We have enforced the RS fine tuning condition, 



\k\1 2 



6, 



(13) 



which means that there is no net cosmological constant 
on the brane. We assume a single component perfect 
fluid for the brane matter, with equation of state p = wp, 
which implies p oc a~ 3 ^ 1+w \ as usual. This allows us to 
write the Friedmann equation as [q = 3(1 + w)]: 



(HI) 2 



(14) 



Here, e* is the energy density of brane matter, normalized 
by the brane tension A, at some reference epoch z b = z*; 
i.e., e* = p*/A. Of course, z* is freely a specifiable length 
scale. 

Now, we turn our attention to tensor perturbations. 
These are defined by the substitution 



1 



J d 3 kh(t,z;k,A)e lkx e 



in the 5-dimensional metric 0. Here, k 



is 



{A) 
ij ' 

(15) 
a 3- 



dimensional wavevector, s 



(A) 



e;r(k) 



is a constant 



transverse trace-free polarization 3-tensor orthogonal to 
k, and the summation is over polarizations. The 5- 
dimensional Planck mass satisfies £k 2 AI 3 = 1. Unless 
explicitly required, we will omit the k and A arguments 
from the Fourier amplitude h below. One finds that h 
obeys 



= 
= 



' dt 2 



d 2 h 3 dh 



k 2 h, 



at oz 



(16a) 
(16b) 



Now, to put these equations into the standard form of 
Sec. Ill Al we just need to make the definition 2 



(17) 



Then, making use of z b 
Eqs. (jTtjj) reduce to 



-H£ and eq. ©, we find that 



= {-D 2 + V(zM 

= [(n-D)i/)-a(ri)il>] b . 
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respectively, with 

V(z) = k 2 + — . 



4z 2 ' 

3 4 _. 3 ^/l + H 2 £ 2 
2 z b 2 z b 



(18a) 
(18b) 



(19a) 
(19b) 



Here, n a is defined by (0. Hence, we have success- 
fully transformed the RS gravitational wave equation into 
an equivalent form defined in the flat (t, z) 2-manifold. 
We call <|18H the 'canonical wave equation' governing RS 
gravitational waves. 

It is convenient to select z* to characterize the epoch 
when a perturbation with a given wavenumber k enters 
the horizon: 



k = a. t H* => kz« = HJ = y/e*(2 + e»). (20) 
If we then work with the dimensionless variables 

i = tjz. t , z = z/z*, fj — r]/z^, k = kz*, (21) 

the brane equations of motion, Eqs. ijHJ and l|14|l . and 
gravitational wave equations, Eq. I|18|) . & re completely 
specified by the single parameter e*. It is useful to define 
a critical value e c = \[2 — 1 « 0.41 that is associated with 
a perturbation that enters the horizon when H*£ — 1. We 
can then classify modes as either long (e* < e c ) or short 
(e* > e c ) wavelength when compared to the character- 
istic length scale £ of the extra dimension. Intuitively, 
we expect the 5-dimensional effects to be important for 
short wavelength modes that enter the horizon when the 
universe is smaller than £. 

Finally, to finishing specifying the problem, we need to 
fix the position of the dil ± boundaries of the computa- 
tional domain in Fig.^ This is equivalent to selecting the 
initial and final brane size: z, = z b (rji) and Zf — z b {rjf). 
It is useful to characterize the initial time by a dimension- 
less number sq, which is the ratio of the perturbation's 



/ 



perturbation enters 
Hubble horizon (a = a*) 




w = 1/3 
= 5 



FIG. 2: (Colour online.) Illustration of how the parameters 
so and cif/a, move the past and future null boundaries of the 
computational domain, respectively. Here, we have taken a 
radiation-dominated brane, so dQ~ is pushed further into the 
past as so is increased. (The opposite is true for a vacuum- 
dominated or deSitter brane, for example.) All other param- 
eters being equal, the amount of CPU time required to com- 
plete one simulation scales with the area of Q divided by the 
typical area of one of the cells shown in Fig. 



wavelength normalized by the horizon size at the initial 
time: 



a,; H.: 



•so 



(22) 



Hence, by choosing so and e*, one determines zi = z^/z*. 
One can place the future boundary of f2 by selecting the 
ratio of the final brane size to the size at the horizon- 
crossing time df/a* — l/if. Therefore, in dimensionless 
coordinates, the computational problem in Sec. Ill Al is 
completely specified by (e*, so, a//a*), up to the choice 
of initial conditions on dVl~ . In Fig. [21 we show how the 
choices of sq and a//a* alter the shape of the compu- 
tational domain for a radiation dominated brane (w — 
1/3). 



III. NUMERIC SCHEME 



2 For our numeric work, it is convenient to characterize GWs by 
the ip wavefunction, as opposed to h. But in the literature it has 
become standard to express perturbations in terms of h, so we 
will always report our results in as h(t, z) instead of ip(t, z). Of 
course, it is trivial to move between the two descriptions using 

im . 



In this section, we develop a numerical scheme for solv- 
ing the class of 1+1 dimensional wave equations intro- 
duced in Sec. Ill Al In subsequent sections, we will apply 
this scheme to the RS tensor perturbation problem de- 
scribed in Sec. Ill Bl 
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A. Computational grid 



composing dO individually. The result is: 



We begin by discretizing the computational domain f2 
into a finite number of 'cells' as shown in Fig.^ Each cell 
is either a diamond whose boundary consists of four null 
segments of equal 'size', or a triangle whose boundary is 
made up of two null and one timelike segment. The cells 
are arranged into rows bounded by surfaces of constant 
u such that all diamonds in a given row have uniform 
size. Also, each row contains one triangle where it inter- 
sects the brane at its leftmost extreme. (In the example 
of Fig. 1, the triangles in the early rows are very narrow 
and hard to see.) The timelike segment in each triangu- 
lar cell gives a straight line approximation to the brane 
trajectory. 

Note that the diamond size generally varies from row 
to row. This is because we have demanded that the u = 
constant row boundaries intersect the brane at evenly- 
spaced intervals of coordinate time t. The magnitude 
of this spacing is given by the parameter A; i.e., the 
future boundary of the i th row intersects the brane at 
ti = to + *A, where to is some initial time. Note that this 
is not the only possible choice; for example, we could 
have demanded each row be regularly spaced in it, 77, or 
some other coordinate. This particular choice of spacing 
ensures that each diamond cell has an area less than 2A 2 , 
while each triangular cell has an area less than A 2 /4. 



B. Diamond cellular evolution 

We now derive formulae that relate the value of the 
field at each node in a particular cell — these will form 
the heart of the integration scheme introduced in the 
next subsection. First, consider a typical diamond cell 



shown in Fig. 3(a) This type of cell has four nodes that 
we label by their compass directions: north, south, east, 
and west. The value of the field at each node is denoted 
by ip n , ip s , xp e , and tpw, respectively. We integrate the 
wave equation Q over the cell and use Gauss's law to 
obtain: 



f (n<> • D) i> = f 

J dO JO 



Vip. 



(23) 



Here, n<> is the outward pointing normal to the cell 
boundary <90. In each integral, the natural 'volume' ele- 
ment on the respective submanifolds is understood. Be- 
cause all of the boundary surfaces are null in a diamond 
cell, n<y is actually everywhere tangent to dO. If A is 
an affine parameter along each segment of the boundary, 
then we obtain 



dx 

~d\ =* (n < 



D ^=dx- 



(24) 



This implies that the boundary term can be evaluated 
exactly by integrating over the four null line segments 



(n ■ D)ip = 2(ip e + ip w ~ip n - ip s ). 



(25) 



dO 



Note that when dealing with each segment, it is impor- 
tant to integrate in the direction of increa sing affine pa- 
rameter A, which is indicated in Fig. 3(a) by the interior 
arrows. 3 

By adopting a bilinear approximation for the inte- 
grand, the volume integral in (|23|) is given by 



-(V n ^ n + V s ij; s + V e ^ e + V w ^ w ) + 0(S 4 ). (26) 



Here, V n is the value of the potential at the northern 
node, V s at the southern node, and so forth; hence, V n = 
V s . Here, 6 is the size of the cell in null coordinates. Our 
choice of grid spacing implies that 

6<2A =s> 6 = 0(A). (27) 

We now use l|23|) and (|2*S|l in (|2*3j> to isolate ipn- 

Tpn = -^s + {lp W + ^e)(l ~ \5 2 V S ) + 0(A 4 ). (28) 

Given the field value at the southern, eastern and western 
nodes, this formula gives us the value at the northern 
node correct to cubic order in A. 

Note that in order to derive l|28|l , we performed a series 
expansion in V s 5 2 and retained the first correction term. 
Hence, one should only believe the C(A 4 ) error term in 
the diamond evolution law when this approximation is 
valid; i.e., when 



8 2 V{z) < A 2 V(z) < 1. 



(29) 



This is sensible requirement: In order to achieve reliable 
results, the characteristic size of a diamond A must be 
much smaller then the characteristic length scale defined 
by the potential 1 / yJV{z). 



C. Timelike triangle cellular evolution 

We now turn our attention to the timelike triangular 
cells at the end of each row, an example of which is shown 
in Fig. 



3(b) This type of cell has three nodes: north, 



south and east, and the boundary is composed of one 
timelike and two null line segments. By construction, 
the brane's trajectory precisely intersects the northern 
and southern nodes. In the calculation that follows, we 



3 Also note that this result could have been obtained via the ex- 
plicit double integral J^D 2 ^ = J J du dv (— 2d u d v ip). This is 
how this formula is usually derived for numeric problems in black 
hole perturbation theory |2fil , but it is difficult to generalize such 
an approach to triangular cells. 



(a) Diamond 



(b) Timelike triangle 



FIG. 3: (Colour online.) Cellular geometries 



view the timelike cell boundary 9Aj as interchangeable 
with the actual brane trajectory between these nodes, 
which we call Of course, there is some degree 

of error associated with such an assumption, so we use a 
' sign to indicate equations that are strictly true only 
dAb. Below, we discuss under which 
signs can be regraded as equali- 
ties. 

Integrating the wave equation over a triangular cell 
and applying Gauss's law yields: 



when <9f2[ s ~~*™' ) = o. 
circumstances these 



(n A • D) $ 



Vi/j. 



(30) 



As before, the null portions of the boundary integral can 
be evaluated exactly. For the timclike segment, we re- 



place the path 9A& with dil 



(s—>n) 



This yields 



/"(nA--D)^«2^ e -V„-^- J{n-D)i>. (31) 

<5A an<~ n > 

Again, we label the field values by the node they are as- 
sociated with. Using the boundary condition JHJ , we can 
substitute for the normal derivative in the brane integral. 
Using the trapezoid approximation, we have 



(n-£))V = ^KVn + « s ^) + 0(^). (32) 



am 



Here, a n and a s are the values of 01(77) at the northern 
and southern node, respectively. The volume integral 



over the cell is handled via a bilinear approximation as 
before: 



A 



Vi>= ^(V n ^ n +V S A+Ve^e)+0[(S 2 U +5 2 V ) 2 }. (33) 



Now, let us denote coordinate differences between the 
northern and southern nodes by St, 8 Zl etc. Our choice 
of grid spacing then gives 



1 1 Si « 8 2 t - 51 < A 2 

2(S 2 + 2S 2 Z )< 6A 2 . (34) 



-jy ■- J t 

5 Z <A 



When we put Eqs. (|3U|I - H34|I together, we obtain 
12 + &a s 8 v + 5 U 8 V V S 



ip-n 



12 + 6 a n 8 n + 8 U S V V 7 

24 - 8 u 5 v V e 
12 + 6 v, 



(35) 



Given field values at the southern and eastern nodes, this 
formula gives ip n accurate to quadratic order in A, pro- 



vided that the discrepancy between <9f2 



( s —>n) 



and dAh is 



negligible. 

Under which circumstances can we ignore this discrep- 
ancy and replace the above signs with '='? Clearly, 



when dfl 



(s->n) 



is well approximated by a straight line 



throughout the cell. In other words, when the change Su 
in u over the cell is small. Note that because the length 
of u is conserved, Su must be parallel to n, so we really 
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want | n • du\ <C 1. This condition can be rewritten as 
|n • <5tx| ~ \n a uP Df3u a \5 ri 
< \u a u D n a \A 
= |D-n|A< 1. (36) 

This condition can be cast in a more geometric light by 
noting that in 2 dimensions, the local radius of curvature 
of a curve is given by r c = r c {rf) = l/\D • n\. Hence, the 
above condition is equivalent to 

A « r c ( V ). (37) 

In other words, we can reliably use the triangle evolution 
law H35fl if the radius of curvature of the brane is much 
larger than the characteristic size of the cell. 

D. The algorithm and theoretical convergence 

Having defined our computational grid in Fig. ^ and 
derived the diamond and triangular evolution laws, l|28l) 
and l|35(l . our algorithm for solving the wave equation 
JTJ is quite straightforward. Referring back to Fig. ^ 
we see that by specifying the value of ip on Q ~, we gain 
knowledge of the field at the past boundary of row 1. 
An enlargement of the situation is shown in Fig. To 
obtain ip on the future boundary of the row, we first use 
the triangle law to obtain the value at the node marked 
'1' in the diagram. Then, we use the diamond evolution 
formula to fill in node '2', then node '3', and so on. 

After the field values on the future half of row 1 are 
found, we then need to find the field at the nodes on 
the past boundary of the next row. As can be seen in 
the diagram, those nodes do not line up with the future 
nodes of row 1, so we must use an interpolation scheme 
to determine the field there. We use a polynomial in- 
terpolation with a four-point stencil: That is, for the i th 
node on the past of row 2, we use the four closest nodes 
on the future of row 1 to find a cubic polynomial approx- 
imation to ip. That approximation is then used to 'fill-in' 
the value of ip at the node on the past boundary of row 
2. This introduces an error of order C(A 4 ) in ip. Once 
this is accomplished for all the past nodes on row 2, we 
then repeat the cycle by using the evolution laws to find 
the field on the future of row 2, then interpolating to get 
ip on the past of row 3, etc. 

Since we have errors of order A in the evolution of 
bulk cells and our interpolation from row to row, and er- 
rors of order A 3 in cells bordering the brane, we expect 
our overall algorithm to exhibit quadratic convergence 
overall, provided that the conditions (|29|) and (|37|l are 
met. That is, if the 'exact' solution of the problem is 
given by inexact while our numerical solution with toler- 
ance A is ipA, we expect: 

ipA(t, z) - V>cxact(i, z) = A 2 e{t, z). (38) 

Here, e(t, z) is a function that does not depend on A. We 
will test this convergence condition explicitly in the next 
section. 




FIG. 4: (Colour online.) The algorithm used to deduce the 
field value on the future boundary of a row, given the value on 
the past boundary. The numbers indicate the order in which 
the nodal fields are calculated, while the arrows inside the 
cells show how information flows through the diagram; i.e., 
the field value at nodes where arrows end depends directly on 
the field value at the nodes where the arrows start 



IV. CODE TESTS 
A. Inertial (de Sitter) branes 

In this subsection, we specialize to RS braneworld cos- 
mological models for which an exact solution to the ten- 
sor perturbation problem of Sec. Ill Bl is known. Namely, 
we consider de Sitter branes with w = — 1. Our goal is to 
compare the results of the numerical scheme introduced 
in the previous subsection to this exact solution, and thus 
test the reliability of our algorithm. 

It is convenient to introduce a new coordinate system 
to describe this scenario: 

£) — ?ycosh£ + Zi coth^, 2(77,^) = — yysinh^. (39) 

Here, £b and Zi are arbitrary positive constants, the time- 
like coordinate rj is strictly negative, and the spacelike 
coordinate satisfies £ > Then, when w = — 1 the 
brane equations of motion @ and (|14|) are solved by the 
£ = £h hypersurface; i.e., 

h{n)=t{riM), z b {r,) = z(ri,t b ). (40) 

The Hubble constant on the brane is given by HI = 
sinh^f, and we find that the brane's speed, 

^ = -tanh6, (41) 
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t = z l coth£. 




£ = const. 



const 



FIG. 5: (Colour online.) The level curves of the Rindler- 
like (77, £) coordinates in the (t,z) plane. A brane undergoing 
pure-de Sitter inflation can be identified with any £ = constant 
surface. The solution to the tensor perturbation problem of 
Sec. Ill Bl is known exactly in these coordinates, which allows 
us to test the accuracy and convergence of the numeric algo- 
rithm developed in Sec. lIIII An example of the computational 
domain Q used for the numeric calculation is also shown. 



is constant; i.e., the brane trajectory in the (t,z) refer- 
ence frame is a straight line. One could easily call such 
branes 'inertial'. Also note that when Zb — Zi, we have 
tb = 0. In Fig. [S] we show the coordinate lines of the (77, £) 
patch in the (t, z) plane. The observant reader will note 
that they are identical to the coordinate lines of Rindlcr 
coordinates in Minkowski space; i.e., the £ = constant 
lines represent a family of inertial observers whose tra- 
jectories are converging to a point and the 77 = constant 
curves represent the surfaces over which their clocks are 
synchronized. 

The advantage of introducing the (77, £) coordinates is 
that they render the wave equation and boundary con- 
dition separable. Hence, it is relatively straightfor- 
ward to obtain an exact solution for h = h(r],£). In 
general, this is given by a superposition of mode func- 
tions {</>o,0„}, with v > 0. The simplest of these is the 
so-called 'zero-mode', 



C(sinh£b) sinh£f, 



,-ikr] 



(42) 



where 



C{x) = 



\/l + x 2 + x 2 In 



1 + VI + x 2 



-1/2 



(43) 



Since h — <fio is a legitimate solution of (|16() . then 



3/2 

Z3/2 



Re 



{[kr ) (t,z)-i]e- ikr >( t < z A (44) 



is a solution of the canonical wave equation (|18|) , with r\ 
defined by 



'1 



(t, z) = -y/(t- Zi COth£fc) 2 - z 2 



(45) 



We also define /i cxact = (z/z*) 3/2 ?/> OX act- 

To test the algorithm, we proceed as follows: We set 
the computational domain Q (shown in Fig. |5J) by fixing 
a de Sitter brane trajectory, and then selecting an initial 
and final time. In practice, we use the dimensionless co- 
ordinates H21fl . so and the dimensionless wavenumber 
fcz* are determined by choosing (e*, so, ay /a*). Then, our 
gird is defined by the selection of a spacing A. We syn- 
chronize our numeric solution i/^a to the exact solution 
on the initial time hypersurface: 



ip A (dtt ) = -0cxact(<9f2 ). 



(46) 



Next, we use algorithm of Sec. IIII Dl to obtain ipA 
throughout fl, and then multiply by z 3 / 2 to get /ia- We 
define the 'distance' between two arbitrary functions on 
the brane as 



1 



1)f - Tji 



(fl - hf 



1/2 



(47) 



which can be thought of as the root-mean-square (RMS) 
deviation between /1 (77) and f% (rj) . Eq. I|38() then predicts 



cr b (A) = ((h A - h, 



exact// b 



const, x A 



provided that 



A < r c and fcA < 1. 



(48) 



(49) 



In the second inequality, we have used that V(z) — 0(k 2 ) 
throughout most of Cl. Note that because de Sitter branes 
have inertial trajectories, r c — 00 and the first inequality 
is trivially satisfied. 

We have calculated Gb for a wide variety of simulations 
of de Sitter brane scenarios and show our results in Fig.|S] 
Each curve represents families of simulations where A is 
varied, but all other parameters are kept fixed. We see 
that (76 is indeed proportional to A 2 for A sufficiently 
small, which allows us to draw two conclusions: First, the 
numeric solution is indeed approaching the exact solution 
in the limit of A — > 0; and second, the rate of convergence 
is quadratic. Notice that this quadratic convergence sets 
in for lower values of A/z» as e* is increased; this is 
because of Eq. l|2Tj|> . which says that fcz» scales like e* 
when e» is large. Hence, in order to satisfy kA <C 1, A 
must approach as e» — > 00. 
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lo gio( A A*) 

FIG. 6: (Colour online.) RMS deviation between exact 
and numeric solutions log 10 ab versus log 10 (A/z*) for several 
simulations of tensor fluctuations around inertial (de Sitter) 
branes. Quadratic convergence, er(,(A) <x A 2 , is evidenced for 
A sufficiently small, indicating our numeric algorithm is both 
stable and accurate under these conditions. 



B. Non-inertial branes 

The code tests in the last section were for manifestly 
non-accelerating branes, so one might reasonably worry 
about the reliability of our algorithm for more compli- 
cated brane trajectories. However, one cannot test the 
numerics in the same manner as before, precisely be- 
cause no convenient exact solution to the wave equation 
is known for an accelerating brane. Indeed, this was one 
of the main motivations of developing the algorithms of 
Sec.InD 

In the absence of an exact solution, we can test for the 
convergence of our numeric results, but not the accuracy. 
In other words, we can confirm that our results stably ap- 
proach some limit (in the Cauchy sense) as A — > 0, but we 
do not know if it is the right answer. A convergence test 
can be formulated as follows: Keeping everything else 
constant, we run our code once with accuracy V2A and 
again with accuracy A (the second run will take around 
twice as much CPU time as the first). Then, we can de- 
fine the RMS discrepancy on the brane between the two 
runs as 

C 6 (A) = ((h A - h A/V2 )) b . (50) 

As in the last section, eq. <|38H predicts that this statistic 
should obey 

Cft(A) = const, x A 2 , A < r c and kA < 1. (51) 

In Fig. we plot ^ versus A for a number of differ- 
ent non-inertial brane trajectories corresponding to radi- 



-2.6 -2.4 -2.2 -2 -1.8 -1.6 -1.4 
lo gio(A/z») 

FIG. 7: (Colour online.) The convergence statistic log 10 
versus log 10 (A/z*) for several simulations of GWs about non- 
inertial (radiation-dominated) branes. As in Fig. E] quadratic 
convergence, 0>(A) oc A 2 , is apparent as A — > 0. 



ation dominated brane universes. The initial condition 
for these simulations is simply 

/i(<9fT) = 1. (52) 

Note that, to a large degree, the convergence properties 
of the algorithm will be independent of the initial data. 
We again see quadratic convergence in this figure for A 
small enough, and that smaller values of e* lead to faster 
convergence. 

From this, we can infer that our algorithm is conver- 
gent when the brane is accelerating. However, our igno- 
rance of the exact solution means that its accuracy is still 
(technically) in question. The only way to get a handle 
on the latter is to compare our results with those inde- 
pendently obtained by some other group/method. This 
is done in Appendix IVII where we explicitly co mpare p ur 
results for a particular brane model to those of HKT. 



V. SETTING THE INITIAL CONDITION FOR 
BULK FLUCTUATIONS FROM INFLATION 

Having convinced ourselves that the numeric algorithm 
developed in Sec. IIIII is returning reliable results, we 
turn attention to the physical problem we want to study. 
This is the the evolution of tensor perturbations in the 
high-energy radiation regime after the end of inflation 
in braneworld cosmology. In this section, we discuss the 
difficulties associated with determining the initial condi- 
tions of our classical problem from a quantum mechanical 
analysis of inflation. 
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(t = +00) 




00) 



transition 
point pt 



FIG. 8: (Colour online.) Conformal diagram illustrating the 
causal structure of a braneworld model of the early universe 



The simplest picture of RS brane cosmology in the 
early universe assumes that the brane has some initial 
phase of pure-de Sitter inflation followed by a period 
of radiation-dominated expansion. These two distinct 
brane trajectories are smoothly joined at some transition 
point pt in the (t, z) plane. The situation is illustrated 
in Fig. [8] via a conformal diagram. To generate this plot, 
we have applied the standard compactification to the the 
(t, z) coordinates: 



T = 



tarda v + tanh u 



tanh v — tanh u 



(53) 



We see that the brane begins in the vicinity of past time- 
like infinity i~ , reaches the transition point at some fi- 
nite (t, z), and then continues expanding as a Friedmann 
brane on to future timelike infinity i + . 

Now, the standard paradigm concerning braneworld 
GW perturbations is they are generated quantum me- 
chanically during the de Sitter phase of the expansion. 
The calculation relies heavily on the (77, £) Rindler-like 
coordinates introduced in Sec. IIV Al since the exact so- 
lution of the classical equation of motion is known ex- 
actly in that patch. The coordinate lines of this patch 



have been drawn on Fig. [SJ One assumes that h(r], £) is 
described by a quantum field in its vacuum state, as mea- 
sured by observers travelling on £ = constant slices. We 
denote this 'de Sitter invariant vacuum' by |0) I( . Then, 
one simply evaluates the vacuum expectation values of 
the squared amplitude of the various mode functions to 
see how they are magnified during inflation; for example, 
the zero mode amplitude is 

,(0|^|0)„ = \MV,0\ 2 = function of V only. (54) 

That is, the evolution of the quantum fluctuations is de- 
termined by the amplitude of the classical mode func- 
tions, which are known exactly. In this way, Langlois 
et al. [l3 have shown that while the amplitude of the 
zero mode 4>q grows during inflation, the amplitudes of 
the so-called 'massive modes' <j} u are suppressed, which 
leads to a scale-invariant spectrum that agrees with the 
4-dimcnsional result, 



1^10), = \Mv)\' 



n C 2 (sinh£ b ) sinh £ b 
~o* 2k 3 £ 3 



Cf(k). (55) 



This quantum calculation suggests that the appropriate 
post-inflationary initial condition for classical GW calcu- 
lations is 



inflation: = const., 



(56) 



where Y,^ is an r\ = constant hypersurface that intersects 
the brane 'at the end of inflation' (cf. Fig. |SJ, which for 
our purposes can be taken as the transition point pr- 

But there is problem with this picture: In order to have 
a well-defined Cauchy problem for the fluctuations h(t, z) 
from the transition time to the infinite future, we must 
have initial data on the entire S„ hypersurface, which is 
the u — constant line running from pt to X + . Now, by 
specifying initial data on we can immediately use the 
bulk wave equation to obtain the field value in it's causal 
future J + (T, V ). 4 But this is not sufficient to determine 
the value of the field on £„; it is apparent from Fig. |S| 
that 



(57) 



Here, an overbar indicates the closure of a set. The rea- 
son that S u does not lie in the future domain of depen- 
dence of H v is the existence of a Cauchy horizon H + in 
the (f?,£) coordinates. In Figs. [S] and OH this is the future 
boundary of the portion of the (t, z) plane covered by 
the de Sitter coordinates. Note that this problem is not 



4 This is possible to do analytically, since the general solution of 
the bulk wave equation is known in closed form, and evolution 
in </+(£,,) proceeds independently of the brane boundary con- 
dition. This holds for any spacelike hypersurface whose causal 
future does not intersect the brane. 
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in general mitigated by choosing to model the brane in 
a finite computational domain O; as in Fig. |H| one can 
draw many examples where dQ~ ^ J + (£,,). 

There are several ways to get around the fact that 
J + (T, V ) is too small. In the literature, several authors 
note that the wave equation and boundary condition Qlfijl 
can be solved approximately for extremely long wave- 
lengths: 

k — ► 0, h(t, z) ~ constant for all (t, z). (58) 

The conclusion that one draws is that, in the limit, the 
inflation initial condition (|56[) is 'consistent' with setting 
h = constant on som e other, more suitable hypersurfacc. 
Consequently, iHKTl have adopted the initial conditions 

IHKTl h(E t ) = const., fo, t (S t ) = 0, (59) 

where E 4 is the hypersurface running from the transition 
point pt to spatial infinity, as shown in Fig. |H1 This is 
sufficient to specify the Cauchy evolution of h, since 

E u C J+(E t ). (60) 

Note that if the brane were static, this initial data would 
reproduce the zero mode of the original RS model: 

fo(t,z) = const, x e~ lkt . (61) 

Of course, this is not a unique prescription; llNl have in- 
stead elected to enforce: 

llNl : /i(E„) = const.; (62) 

i.e., they have set the perturbation equal to a constant on 
the initial null hypersurface. Both groups acknowledge 
that these initial condition are somewhat ad hoc; when 
k = they are only consistent with the inflation initial 
condition Ij56(l if additional data is specified on T~ , and 
they are not even consistent with one an oth er fo r k > 
0. We will explicitly contrast the IHKTl and IN initial 
conditions in Sec. I VII 

A separate approach comes from treating the quan- 
tum inflationary calculation differently. Gorbunov et al. 
[20| consider a 'junction model' that has de Sitter and 
Minkowski branes attached to one another at a non- 
smooth transition point. They assume that the GWs 
are in the de Sitter invariant vacuum |0)^ in the infinite 
past, which implies quantum particle creation when the 
brane's trajectory changes abruptly. At the end of the 
day, then derive the spectrum of GWs as seen by RS ob- 
servers travelling orthogonal to t = constant slices. At 
first glance, this would seem to be ideal since the final am- 
plitudes are given on E t ; furthermore, the actual derived 
spectrum is domin ated b y the zero mode, and is hence 
consistent with the lHKTl initial condition (|59|) . However, 
some caution is required here, since the calculation essen- 
tially involves decomposing the quantum states preferred 
by z = constant observers in terms of those preferred by 
£ = constant observers. But the latter basis is only really 



defined to the past of the Cauchy horizon, so we may le- 
gitimately worry if the field amplitude on the entirety of 
Ti t is fixed in this approach. (This issue is clearly beyond 
the scope of the current work.) 

Recently, Kobaya shi and Tanaka [lj| have modified 
the iGorbunov et alJ calculation by smoothly joining the 
initial and final phases by a radiation-dominated brane. 
In order to calculate the Bogliobov coefficients governing 
particle creation, they follow the usual practice of fixing 
the field value in the future (on an v = constant slice) and 
calculation what it is in the past (on an initial v ~ con- 
stant slice like E„ in Fig. |SJ). They found a subdominant 
amount of energy < 10% is stored in the Kaluza-Klein 
modes at the end of the quantum regime. Kobayashi 0] 
subsequently repeated the calculations with an improved 
numerical scheme, which resulted in more robust predic- 
tions. 

Finally, we note one additional strategy for specifying 
initial data; namely, to enforce boundary conditions on 
T~ . Imposition such an additional constraint is enough 
to fix the problems associated with the inflation initial 
condition (|56|l . since 

E„ C j + (E,ur). (63) 

The real question is: what is the appropriate condition to 
choose? A physically sensible minimal condition would 
be that there is no radiative flux entering our patch of 
AdS space through Z~ . 5 It would be interesting to exam- 
ine this type of prescription in detail, but such a project 
is beyond the scope of the current work. 



VI. HKT VERSUS IN FORMULATIONS 

In the previous section, we saw several different types 
of initial conditions that one could use for numerical sim- 
ulations of bulk GWs in the post-inflationary epoch of 
brane cosmology. In this sectio n, we e xplicitly compare 
the GWs generated by the HKT and IN initial conditions. 
It is useful to first identify the key qualitative feat ures of 
the radiation, and to do this we concentrate on the[HKT 



5 One could argue that, in some sense, the potential enforces this 
condition for us: Far from the brane, if) behaves like a field of 
mass k propagating in free space. Hence, wave packets of if) 
must originate at i~ , and not from past null infinity. In other 
words, wave packets on I~ should never reach the brane, and 
therefore it seems that initial data on T~ can have no relevance 
to the observed GW spectrum. However, this conclusion relies 
both on the properties of the potential (and hence fails for k = 0), 
and an eikonal approximation. Here, we are interested in making 
statements independent of the potential and any approximations, 
which means that specification of initial data on X~ is logically 
distinct from the specification of data on E^, E„, etc. Also, 
from a practical point of view, it is of little use to assume that 
our solution is independent of the field configuration on X - ; our 
numeric code crashes unless data is specified on an initial u-slice, 
no matter how far in the past it is. 
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case. For simplicity, we take the HKT condition to be 
h = 1 and d t h = on an initial t slice in actual calcu- 
lations. Labeling the initial time as t = 0, we can write 
down the exact solution of the wave equation l|l(j|) inside 
the causal future of S t : 

h(t, z) = cos(fci) for all (t, z) £ J + (£ t ). (64) 

In particular, this can be evaluated on to give us 
initial data on <9f2~, 

h(dfl') = cas±k(v + Ui), (65) 

where u% is the u-coordinate of the intersection of the 
brane with the initial time slice. 

In Fig. ED we plot the results of our numeric calcula- 
tion for a particular case involving a radiation dominated 
brane. The key features are as follows: In the bulk, the 
waveform appears to retain the character of the Randall- 
Sundrum zero mode far away from the brane; i.e., it is 
constant on constant t slices. However, this symmetry is 
broken closer to the brane by the motion of the bound- 
ary, resulting in rich and intriguing dynamics. On the 
brane, the GW amplitude remains constant for a < a*, 
and then begins to oscillate and decay. Asymptotically, 
one can easily confirm that the numeric waveform goes 
like 

a C b (k) ( a \ C b (k) 

rib ► cos I Wo h? = cos(fc?7 + <;), (66) 

oo a \ a* / a 

where Cb(k) > is the expansion-normalized character- 
istic amplitude, <; is a phase, and 

lo = y/l + %Z. (67) 

This asymptotic form can be understood by considering a 
fictional 4-dimensional universe that undergoes the same 
expansion (c/. Ea. ll4|) as our 'real' brane universe. Tensor 
fluctuations in such a model obey 

( + 2Ha4- + k 2 ) h ief = 0, (68) 
\ drj z dr\ J 

where the 'ref subscript indicates the amplitude in our 
reference 4-dimensional universe. At late times, Hi ~ 
V / 2e*"(a*/a) 2 for a radiation dominated brane, leading to 



HKT formulation 



h(E t ) = lA(£ t ) = 

t . 4-dimcnsional 

late time behaviour 




(a) Bulk profile 



(e„a ) = (5,10) 




-0.2 0.2 0.4 0.6 0.8 1 1.2 
log 10 (a/o*) 



(b) Brane profile 



h Te {(a) 



cos(fc?7 + w), C rcf (fc)>0. (69) 



Hence, the late time behaviour of h b from the numeric 
simulations matches that of the 4-dimensional effective 
calculation. However, the amplitude of h b relative to h rc { 
is somewhat diminished, as seen Fig. 9(b) Physically, 
this has to do with the brane motion inducing the radia- 
tive loss of GW energy into the bulk, as depicted by the 



bulk gravity wave profile seen in Fig. 9(a) 



In Fig. 1101 we plot o ur h b predictions for several ad- 
ditional cases using the lHKTl initial condition. For com- 
parison, we also plot the actual numeric results obtained 



FIG. 9: (Colo ur online.) Results of a typical integration us- 
ing the H KT] initial condition (A/*, = 2~ 10 ~ 10~ 3 ). In 
Fig. |9(b)| we have drawn what the brane GW signal would be 
if the bulk were neglected; i.e., if one solved the 4-dimensional 
master equation 1681 with a modified expansion rate given by 
(1121) . Note that we have enforced the initial condition that 
ht = h Ic f for a < a, on the reference wave. 
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FIG. 10: (Colour online .) Direct comparison of numerical re- 
sults obtained bv lHKTl and the current work (A/z* = 2~ 10 ~ 
10 -3 ). In each panel, the RMS deviation between the two 
profiles is given. 



bv iHKlj for the same cases. The agreement between the 
two independent calculations is visually excellent, with 
both sets of data lining up very well, but not perfectly. 
We can quantify the level of agreement by calculating the 
RMS discrepancies ((/ia ~ ^hkt))6, which are written di- 
rectly on the figures. We see that ({Ha — ^hkt))& ^ 10~ 2 
for the cases shown, which can be interpreted as the av- 
erage absolute deviation between the two results. This 
is acceptably small, and we can conclude that the two 
simulations agree to within reasonable tolerances. 

Now, we turn our attention to the lINI formulation. In 
Fig. 1111 we show the radiation patter ns arou nd a w = 1 
'stiff matter' brane generated by the iHKTl and |S| ini- 
tial conditions, respectively. The two bulk waveforms in 
Fig. 11(a) are similar to one another, but exhibit some 
clear differences, especially near the initial data s urface 
dfl~. But the the brane profiles shown in Fig. 11(b) 
are virtually identical to one another. Indeed, the RMS 
discrepancy between the two results is ~ 1CP 4 , which is 
quite small. 

Is it possible t hat th e re mar kable on-brane agree- 
ment between the HKT and [S| formulations is due to 
a serendipitous choice of parameters? In Fig. ^1 we plot 



HKT formulation IN formulation 

t /N /i(Ei) = l,/i t (E t ) = u h(E„) = l 




(a) Bulk profile 



o 
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offset 
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(e„s ) = (1,3200) 
<</4 HKT) -fc A N) >>6 = 6-6x10 
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(b) Brane profile 



FIG. 11: (Colour online.) GW a mplitud es h A (t , z) for a 'stiff 
matter' brane (w = 1) using the IrlKTl and ITnI initial condi- 
tions. 



the discrepancy (as defined by the inner product Eq. I47|) 
between the brane amplitude generated by the two dif- 
ferent initial conditions as a function of the initial wave- 
length of the perturbation sq- In this plot, the brane is 
radiation-dominated. We see the limiting behaviour 



«4 HKT) 



(70) 



i.e., as the initial data surface dfl~ is moved further and 
further into the past (cf. Fig.|5J), the level of agreement 
between the two boundary condit ions in creases. This is 
reminiscent of a result obtained by H Klt Their on-brane 
waveforms showed very little dependence on sq , provided 



14 



S" -1 



tuO 

,2 -3 



a 

o 



-4 




0.5 



1.5 



2.5 



Initial wavelength log 10 sq 



FIG. 12: (Colour onlin e.) RM S disc repancy in the hi, wave- 
forms generated by the lHKTI and llNI initial conditions versus 
the initial wavelength log 10 so of the perturbation. 



that So was large enough. 

As described in Appendix 1X1 we need to know 



K 



Cb_ 

C Te f 



late time amplitude of hb 
late time amplitude of h 



rcf 



(71) 



in order to predict the observable spectral density of GWs 
today. Here, we will take the post-inflationary epoch 
to be purely radiation-dominated. By 'late time', we 
mean that the ratio should be evaluated in the low energy 
regime of the cosmological evolution: 



1 > 



(t)' 



(72) 



For practical calculations, we measure the amplitude ra- 
tio when a = 20a* and limit e* < 32. In Fig. 1 1 3 | we 
plot 7Z as a function of f/f c for the iHKll and llNl ini- 
tial conditions and several choices of so; where / is the 
present-day frequency of the simulated wave, and f c is 
the present-day frequency of a mode that re-entered the 
horizon when H£ — 1 (cf. Ea. lA12l) . Qualitatively, we see 
that for so 1, there is a discernable difference in the 
7Z predictions from the two initial conditions. However 
for so 2> 1, the ratios are identical to one another, and 
satisfy 



^«(/// c r 2/3 , /»/c 



(73) 



When this result is combined with (|A9I) and (|A14|) . we 
find that the present day spectral density of gravitational 
radiation obeys 



^cw \ ( j- 



f » fc 



(74) 
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FIG. 13: (Colour online.) The ratio 1Z of the late-time ampli- 
tude of GWs on the brane obtained from our simulations to 
that of the reference wave /i ro f. In all cases, we have assumed 
radiation domination w = 1/3 and evaluated the ratio when 
a = 20a, , or equivalently p/X < 2 x 10~ 4 . 



i.e., we get a flat GW spectrum in the high-frequency 
limit. This essentially re-confirms the main result of 
IHKTl but with the twist that it also holds for the ini- 
tial data prescription favoured bv llNt However, the lat- 
ter group claims the spectrum is red at high frequencies: 
J]qw tx (/// c ) -0 ' 46 - The source of tension between this 
and the current result is unclear. 



VII. GENERIC INITIAL DATA 



In the previous subsection we saw that if our initial 
data surface was set far enough into the past, the on- 
brane wa vefo rms were insensitive to the choice between 
HKT and llNl initial conditions. However, it is clear that 
either alternative represents a fairly restrictive choice of 
initial data. In this subsection, we will explore the extent 
to which hb is indifferent to more arbitrary choices of 
h(d£l~~ ). Throughout this entire section, we specialize to 
radiation dominated w = 1/3 models. 
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A. Basis functions 

Generically, any initial data on <9f2~ may be decom- 
posed as 

oo 

h(m~) = hi{v) = i- J dfiA(fi)e^ v+ ^/ 2 , (75) 

— oo 

where the Fourie r amplitude A(/j,) is arbitrary. For ex- 
ample, the lHKll condition is recovered if 

A(ji) = ir[5(ji-l)+5(n + l)], (76) 

while the |Sj condition follows from 

AQi) = 27t5( M ). (77) 

Let us now define Xm(^) ^° ^ e a solution of the wave 
equation such that Xn{d^~) = e ^k{v+m)/2_ Then; 

OO 

M^) = ^ / i/i^(/x)xi.(n); (78) 

— oo 

i.e., if we have knowledge of the 'basis functions' Xp,i we 
can write down the solution for h corresponding to ar- 
bitrary initial data. Our present goal is to numerically 
calculate x M for various situations and gain some intu- 
ition about its qualitative behaviour. 

Before proceeding, it is useful to make a few remarks 
about {x/j}- First and foremost, this basis is in no way 
preferred or special, it is merely convenient. Its definition 
is intimately tied to the choice of the initial null hyper- 
surface; hence, if c?fl is moved, we get a different basis. 
We do not expect {x/x} to satisfy any type of orthogonal- 
ity relationship over fi; indeed, we will not even specify 
what the appropriate inner product might be. 

An important property of {x^} is the physical inter- 
pretation of the real and imaginary parts of the basis 
functions. If pi is the point where the brane cT2f, and 
the initial data hypersurface <9f2~ intersect, the real and 
imaginary parts of Xy. satisfy 

Re Xtl (Pi) = l, Jmx„(Pi) = 0. (79) 

Hence, the two independent components of any given Xn 
represent distinct physical possibilities: Either the initial 
brane amplitude is nonzero or not. From a 4-dimensional 
point of view, Re Xn represents a superhorizon perturba- 
tion whose non-zero amplitude is frozen until re-entry. 
On the other hand, Imx^ is a perturbation that exists 
'entirely in the bulk' initially; its appearance on the brane 
after the end of inflation would be mysterious to an ob- 
server unaware of the extra dimension. In the special 
case of /i = 1, it is easy to confirm that the imaginary 
part of Xn satisfies the following initial conditions on the 
initial t = constant surface: 

Imx^=i(S t ) = 0, 9Jmx M =i(S t ) £ 0. (80) 



When this is compared to l|59[l . we see that Im^=i sat- 
isfies a sort of 'complimentary' HKT condition. On the 
other hand, Re^=i satisfies the HKT condition pre- 
cisely. 

Finally, we mention the physical interpretation of the 
fi parameter. If we neglect the brane boundary condition 
and work in the limit of z — > oo, we find that 

Xfj, ~ e ik(u-Ui)/2fj, e ip,k(v+Ui)/2 ^g-^ 

is a solution of i|16a[l that satisfies the appropriate ini- 
tial conditions on <9f2 - ; i.e., u — m. It is instructive to 
calculate the flux associated with this 'solution': 

3fi ■ dx = — (x»Dx M - X^DxD ■ dx 

w k(cosh /3dt + smh{3dz), (82) 
where we have defined 

u 2 - 1 

tanh/3=^ . (83) 

pr + 1 

Hence, asymptotically far from the brane and the origin 
(z = 0), the x^ basis functions reduce to plane waves 
traveling with Lorentz boost parameter /3 with respect 
to the static (t, z) coordinates. Note that the definition 
of (3 implies that modes with \p\ < 1 have wavevectors 
pointing towards z = 0, while the m odes w ith > 1 
are traveling away from z = 0. The IHKTl modes /i = 
±1 represent the middle ground: they have no relative 
motion with respect to the static frame. Indeed, when 
(i = ±1 we have Xm ~ e ±lkt in the asymptotic region, 
which are the two independent phases of the RS zero- 
mode. 

To summarize, in this subsection we have introduced 
a set of basis functions in terms of which any square- 
integrable h(dti~) can be decomposed. It must be 
stressed that while this choice is convenient, it is arbi- 
trary. Obviously, if we employed any other basis, the 
interpretation of /i would be quite different; for example, 
we could have selected the RS massive mode functions 
(evaluated on dfl~) from the static brane case as a 
basis, which is perhaps more suited to the fact that the 
bulk is warped. Having said this, our choice of Xn is ex_ 
tremely straightforward to implement numerically, and 
we feel that the /i parameter has an 'easy' physical inter- 
pretation: both as the relative frequency of initial data, 
and simply related to the Lorentz boost parameter of a 
plane wave irradiating the brane from z = oo. 

B. Numerical results 

In this subsection, we present our numerical results 
concerning the evolution of the X n basis functions intro- 
duced above. Note that by introducing this basis, we 
have increased the parameter space from (e*, sq, a//a*) 
to (e*, so, a//a*, /i). Each of these has a continuous spec- 
trum, so it is very impractical to sample this parame- 
ter space densely. Instead, we will attempt to identify 
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FIG. 14: (Colour online.) Bulk GW profiles for IniXn initial 
data. Here, we have taken (e*,so) = (10,5) and w = 1/3. 
The white arrows indicate the direction of the asymptotic 
flux derived in Eq. 



the principal trends in the waveform behaviour from the 
simulations, and go on to rationalize them with some ap- 
proximate analytical work in the next subsection. 

In Fig. ^] we show the bulk GW profile associated 
with the imaginary part of x^ f° r several different choices 
of [i. In all cases, we see that the GW profile far from the 
brane appears to be that of a plane wave, in agreement 
with the discussion of the last subsection. A white arrow 
indicating the expected direction of propagation of these 
plane waves (cf. I82f> is superimposed on each plot; and 
we see that there is reasonable agreement between our 
Xfj, approximation l|81[) and actual simulations asymptot- 
ically. Intriguingly, we see that as /x is increased, more 
of the initial data seems to 'reach' the brane. Stated in 
another way: When fj, is small, the initial data loses its 
coherence when propagating from <957~ to dflb- 

In Fig. 1151 we plot the on-brane waveforms for cases 
similar to those shown in Fig. 1141 Here, we see that Im x M 
is smaller than Rex/i by several orders of magnitude, 
which implies that 



X^(0 6 ) » Rex M (0fi 6 ). (84) 

That is, at first approximation, x M (9r2&) is independent 
of Im^. Furthermore, the real part of x^(d^b) shows 
virtually no variation as fj, is increased. This leads us 
to hypothesize that the brane waveforms are principally 
determined by value of the initial data on the brane: 
The very weak dependence of x^(^^fc) on I m XM i m P ncs 
that the brane signal is insensitive to data with no ini- 
tial amplitude on the brane, while the insensitivity of 
Re% M ((9r2fc) to /i implies that the brane signal does not 
overtly care about the details of initial profile in the bulk. 
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FIG. 15: (Colour online.) On-brane waveforms for Xn initial 
data as /j, is increased and other parameters are kept constant. 
One can see that Re Xm i s relatively insensitive to /i, while the 
opposite is true for ImX/j- Generally speaking, the amplitude 
of the imaginary part is less than the real part by several 
orders of magnitude. 



We now want to describe how the late time amplitude 
of the real and imaginary parts of Xn depend on s and [i. 
But there is small problem: It is apparent from Fig. 1151 
that the late time brane behaviour of Im x^i is much more 
complicated than that of Re Xn ■ Hence, it is more diffi- 
cult to obtain the asymptotic amplitude of the imaginary 
parts without running simulations for a very long time, 
which is computationally expensive. However, if we are 
just interested in a rough characterization of the late time 
amplitude, we can define the expansion-normalized aver- 
age power as 



(a 2 Re 2 X/j)b 



1 



a 2 — a% 



da 



, (85) 



with a similar expression for Im^. For our calculations, 
we select the lower integration limit to correspond to an 
epoch well into the low-energy regime: 



(86) 



On the other hand, is selected to be as large as is com- 
putationally feasible. Roughly speaking, y/2 x (power) 
gives the characteristic amplitudes seen in Fig. 
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Using the average power statistic as denned above, we 
study the dependence of the late time waveforms on pa- 
rameters in Fig. El In both panels, we see that Re \fi 
does not show much variation with either /x or sq for sev- 
eral different values of e*. However, the imaginary parts 
of the basis functions change by several orders of magni- 
tude as the parameters are varied. The principal trends 
are for (a 2 Im 2 x^) to increase with fj, and decrease with 

The physical inferences we can draw from the numeric 
work of this subsections is as follows: The late time be- 
havior of perturbations on the brane is largely fixed by 
the value of initial data at the brane for a large range 
of parameters and choices of hi (v) . The degree to which 
hi, is oblivious to the bulk initial data profile increases 
as the initial data surface is pushed further and further 
into the past (cf. Fig. 01 . Conversely, we find that these 
conclusions are mitigated if the Fourier spectrum of the 
initial data along <951~ involves high frequencies. In such 
cases, the bulk initial data finds it 'easier' to reach the 
brane directly. 

The main observational inference of our simulations 
is that, for \i less than some cutoff fx c , the late time be- 
haviour of the basis functions independent of fi and given 
by 

Xu(dttb) — cosik-q + s), for all ju < \x c . (87) 

oo a 

In general y, c will be a function of So , and we expect that 
fi c — » oo for so — * oo. The exact definition of /i c will 
depend on how rigorously one wants to enforce the ap- 
proximate asymptotic form JS7J). Now, if the amplitude 
components of a particular initial data profile satisfy 

A(n) w for fi > fi c , (88) 

the late time waveform follows directly from (|78(l : 

h b -^-> hAvi) x — cos(fc?y + ?), (89) 
oo a 

where hi(vi) is just the initial data at the brane position. 
By definition, the reference wave introduced in Sec. IVII 
is also linear in hi(vi), so the 5D/4D amplitude ratio TZ 
— and hence f2cw — will be independent of the details 
of the initial data profile, provided that the condition 
(|88|l is met. In other words, for finite so we expect the 
late time gravitational wave spectrum to be independent 
of the detailed shape of the initial data profile, provided 
that that profile does not involve high frequency features. 
From this, it follows that the p red icti on of a flat GW 
spectrum derived from the lHKTi and llNl formulations will 
generalize to more generic initial data. 
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FIG. 16: (Colour online.) Expansion-normalized average 
power of the real and imaginary parts of Xm as functions of pa- 
rameters. To calculate the average, we have selected S = 1CP 3 
and ci2 = 50a,. The average power at late times is not very 
sensitive to these choices. 



C. Analytic results 



The inferences of the last subsection are just that: in- 
formed intuition about the behaviour of initial data based 
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p = (u»,v p ) 



horizon crossing 
p* = («*, v*) 



v(n) « fc 2 




FIG. 17: (Colour online.) Spacetime geometry used to ra- 
tionalize the behaviour seen in the numeric simulations of 
Sec. I Vll^l using an approximate retarded Green's function. 



on the experience gleaned from a finite number of simu- 
lations. We can put them on somewhat surer footing by 
doing some approximate analytic calculations, which is 
the purpose of this subsection. 

Let us analyze the behaviour of the wavefunction ip in 
the high-energy epoch of the cosmological evolution. To 
be more precise, we limit our attention to modes with 
e* » e c w 0.41; i.e., modes that have H r £ » 1. This 
implies that the slope of the bulk brane trajectory when 
a = a* is 



dt b 



(90) 



Hence, the brane trajectory is very nearly null for a < a*; 
here, we will boldly assume that it is exactly null. The 
computational domain f2 for this situation is illustrated 
in Fig. El We label the brane's position at horizon re- 
entry as p* = (it* , u*) , and define II to be the portion of f2 
located to the past of u = u„. Our goal is to approximate 
the field on the future boundary dli + of II given initial 
data on d$l~ . 

Note that since we are assuming that e* is large, then 
fcz* = ^/ e *(2 + e*) is also large; i.e., k > This 
means that the potential in the bulk wave equation 118|) 
satisfies 



V(z) 
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for z > z* 



(91) 



With this, we can immediately write down an approxi- 
mate retarded Green's function for the bulk wave equa- 
tion with (p, p') e II: 



1 



G(p;p') = ^Jo(kX)9(u - u')9{v - v'), 



(-D 2 + k 2 )G(p;p') = 5(p~p'). (92) 



Here, we have defined A = y/(u — u')(v — v') which is the 
proper time interval between the source point p' and the 
field point p. In terms of this Green's function, we can 
express the value of ip at any p £ H as: 

V>(P)« / n.[G(p;p')D'^(p')~^(p')D'G(p;p')}. 

an 

(93) 

Here, the integration is over p', n is the outward point- 
ing normal, and D' indicates the derivative with respect 
to p'; i.e., differentiation with respect to the primed co- 
ordinates. The approximation sign comes from the fact 
that G is not the 'true' Green's function. 

Let us now push p to the future boundary dU + of H 
Since the Green's function has support when p is in the 
future of p', to calculate the integral (|93|l we need only 
specify the field values on <9f2~ and the portion of the 
brane dilb to the past of dH + . We leave the initial data 
on c?f7 arbitrary: 



rz \3/2 
^(dfl-) = (^) hi(v). 



(94) 



However when the brane trajectory is null, the boundary 
condition (I16bll reduces to 



dh 

dv 







3/2 



hi(vi); 



(95) 



i.e., h is constant on the brane before horizon crossing, 
which is entirely consistent with our numeric simulations. 
Using this boundary data, simplifying the integral Q93|) 
is straightforward, but tedious. The result is: 



V>(p) « (2z,) 3 / 2 



3hi(vi) / du 



, J (h\i) 



dv'J (k\ 2 ) 



d hjjv') 
dv' (v' — Ui) 3 / 2 



Oih'JuJ" 1 ) 



-Jo(k\ 



u'f/ 2 

. hi(vj) 

(Vi - Ui) 3 / 2 



(96) 



where we have defined 





- u)(v p - 


■ v*), 




- Ui)(v p 


-v), 



(97) 
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To arrive at this expression, we have used integration by 
parts to remove derivatives of the Green's function. 

The approximation (|96(1 can be used to justify the 
trends we have seen above, if we note that both the ini- 
tial data hypersurface and horizon crossing epoch are in 
the high-energy regime: 



■so 



Hidi 



e* 



O(-ul). (98) 



When inserted in l|96[l . this yields 



+ I dvO{ti iS y 2 ) + 0{hi S y 2 ). (99) 



If we then hold the initial hi(v) profile constant, as in 
Fig. 16(b)| we see that the last two terms drop out as 
so — > oo. Hence, we have shown that as the initial data 
surface goes to the infinite past, the wavefunction on <9II + 
depends only on the value of the initial data on the brane 
hi{vi) — and not on the initial data in the bulk hi(v > 
Vi). Since the evolution of the GWs to the future of dU + 
depends only on the field value on <9IT + , we can conclude 
that the entire late time brane waveform is determined by 
hi(vi) in the limit sq — ^ oo. This explains the behaviour 
seen in Fig. 



Now, consider the situation if we hold the position of 
the initial data hypersurface constant and vary the initial 
data profile, as in Fig. 16(a) Let us also assume that so 
is large and that hi(v) — 0(1). Then, the third term on 
the righthand side of l|99|l is not important relative to the 
first. However, the second term can be arbitrarily large 
if we allow for initial data with large gradients h^(v). In 
particular if we have hi(v) — e *i- lk ( v + u i)/ 2 ^ which is the 
initial data that generates x M , we see that the second 
term grows with /i. Hence, we have shown that if the 
initial data on f?f7 involves a significant high frequency 
component, then the late time brane waveform will be 
sensitive to the precise form of hi(v). This rationalizes 



the behaviour seen in Fig. 16 (a 
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FIG. 18: (Colour online.) Wavelength of perturbations at the 
end of inflation as a function of the inflationary energy scale. 



VIII. CONCLUSIONS 



one-brane scenario. One can find at least three different 
prescriptions in the literature 0, 0, 0] for how to spec- 
ify post-inflationary bulk initial conditions in such mod- 
els, which yield at least two contradictory predictions for 
the spectral tilt of the stochastic GW background gen- 
erated by brane inflation. Until now, it was unclear if 
the discrepancy between various results was due to dif- 
ferent choices of initial data or the numeric scheme used 
to evolve GWs through the high-energy radiation era af- 
ter inflation. Using our code , we h ave inve stigated the 
initial conditions proposed bv iHKTl and lINl and we find 
no discrepancy between the late time GW spectrum gen- 
erated by either prescription, provided that the initial 
data hypersurface is set far enough into the past. Obser- 
vationally, this means that both initial conditions lead 
to a flat GW spectrum flow ^ (///c)° at frequencies 
higher that a threshold f c set by the curvature scale of 
the bulk. This is also in agreement with the results of 
iKTl who study the quantum (i.e. not classical) evolution 
of the GW wavefunction. 



In this paper, we have developed a new numeric algo- 
rithm (from first principles) to deal with the problem of 
solving (l+l)-dimensional wave equations in the presence 
of a moving boundary. Our technique, which is based on 
characteristic integration methods from black hole per- 
turbation theory, has demonstrated itself to be both re- 
liable and accurate by reproducing previously known an- 
alytic and numeric results; and it is the principle result 
of this work. 

We have applied our formalism to the cosmological 
problem of the evolution of GWs in the Randall-Sundrum 



We have also considered more general initial data by 
introducing a practical basis in which to decompose gen- 
eral solutions of the wave equation. Numerically, we find 
that the late time brane waveform is not very sensitive 
to the detailed initial data profile if we start our simu- 
lations at sufficiently high energies. However, this ap- 
proximation breaks down if the Fourier transform of the 
initial data involves very high frequencies. We have used 
an approximate Green's function analysis to analytically 
rationalize these results and demonstrate how they ap- 
ply to any initial data; not just the choices we modelled 
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explicitly. 6 

So, is the high-frequency stochastic GW spectrum pre- 
dicted from this class of braneworld models really flat? 
The answer seems to depend rather intimately on the en- 
ergy scale of inflation. Throughout this paper, we have 
treated so as a free parameter, and our results on the 
insensitivity to initial data depend on the limit Sq — > oo. 
But really, we should fix the energy scale of inflation 
E m { = Hi n {£ and synchronize <9f2~ with the beginning 
of the high-energy radiation era. This means that s is 
actually a function of e* and Ei n f. 



E: 



3/4 



\A*( 2 + e *) 2 



E„ 



1 1/4 



(100) 

For reference, we have plotted Sq as a function of E- m { 
in Fig. E| As can be seen from this plot, if we con- 
sider moderate inflationary energy scales E- m f < 10 3 , it 
is possible to have sq < 10 2 for reasonable values of e*. 
The results of Sec. I VIII imply that such values of so im- 
ply the dependence of the late time waveforms on initial 
data is weak, but non-trivial. This suggests to us that 
the greatest chance of obtaining deviations from the flat 
spectrum lies in models with small inflationary energy 
scales. In such scenarios, it is possible for the brane sig- 
nal to carry a signature of h(dtt~), which in turn de- 
pends on the details of the inflationary model. On the 
other hand, if the inflationary energy scale is high, the 
brane signal will only depend on the value of the pertur- 
bation on the brane at the end of inflation. Investigating 
the detailed gravitational wave spectrum generated by 
moderately low-energy brane inflation is an interesting 
avenue for future work. 
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APPENDIX A: CHARACTERIZING THE 
GRAVITY WAVE BACKGROUND 

In this appendix, we review how to map the results of 
our calculations into observable predictions concerning 
the spectrum of the cosmological GW background today. 
The treatment is very much the same as in Refs. [9t Il2|. 

On the brane, the complete GW perturbation is writ- 
ten as 



hij(r, x) 



(2^Af 5 ) 3 



J2 J d 3 kh b (r;k,A)e ik ^f (k), 



A=+,x 



(Al) 

Here, r is the cosmic time; i.e., ds 2 = —dr 2 + a 2 c?x 2 
on the brane. The energy density associated with this 
perturbation is 



1 / dhfj dh lj 



dr dr 



(A2) 



where the angular brackets indicate an average over some 
region of spacetime. To preform the spatial average, we 
make the standard assumption that the background is 
unpolarized, stationary, and isotropic. This means we 
can trade the operation of spatial averaging over x for 
ensemble averaging at x = 0, and the Fourier amplitudes 
obey: 

(h b (T;k,A)hl(r;k',A')) = 

(2nM 5 ) 3 5 AA ,6(k - k')\h b ( T ; k, A)\ 2 . (A3) 

The product of polarization tensors can then be reduced 



4 A) (k) £ (^(k) 



2. 



(A4) 



Temporal averaging in the late universe can be done by 
noting that in the late universe, our numeric results give 
that hb is a superposition of the e ±lkr> /a 'zero- mode' 
functions. Neglecting scale factor derivatives and using 
r\ « r / 'a yields 



dh b {r;k,A) 



dr 



k 2 C 2 b (k) 
2a 4 : 



(A5) 



where C b (k) is the expansion- normalized characteristic 
amplitude (cf. Eg. I69fl . To relate C b (k) to the primordial 
GW spectrum, we make use of the reference wave /i re f 
discussed in Sec. IVII Our numeric simulations can be 
used to find 



6 Note that this is in general agreement with the recent (numerical) 
work of Kobayashi Il5l using a quantum formalism. But note 
that our result differs in that we co nsider completely arbitrary 
initial conditions, while Kobav ashU considers initial conditions 
from pure de Sitter inflation. 



K = 



C 



rcf 



(A6) 



i.e., the ratio of the characteristic amplitudes in the low- 
energy regime. This is useful because the evolution of the 
reference wave through the high-energy radiation epoch 
is extremely simple: Essentially, it remains constant until 
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a = a*, then its amplitude decays as 1/a. Hence, we can 
write 



CUk) « alCf(k) = (k/H*) 2 C 2 (k). 



(A7) 



Here, C 2 (k) is the squared amplitude of hb, set after infla- 
tion. (Recall that, by definition, hb and h Te { are identical 
before horizon re-entry, which means the share the same 
initial power spectrum.) We can conveniently re-express 
this in terms of 



2 2 



(2nM 5 ) 



(A8) 



which for the inflationary primordial spectrum discussed 
in Sec. W\ reduces to 



6$ = 2 K iC A (H in{ £) 
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(A9) 



where H m { is the Hubble parameter during inflation, and 
we have made u 
together yields 



we have made use of M^Kgi = 1. Putting all these results 



Pgw 



8k 2 4 cl a 



dlnfcfc 2 



k 



s 2 n 2 . 



(A10) 



For comparison to actual experiments, it is convenient 
to re-express this in terms of the frequency / = k/2ira 
observed today and introduce the spectral density 



flGW = 



1 dp, 



GW 



Pc d\nf 



2^f5 2 (f)TZ 2 (f) 
3H 2 H 2 (f) • 



(All) 



where p c = 3H 2 /k 2 is the critical density. 

To progress further, we need to know H* = i/*(/); 
that is, for a given mode with frequency /, we need to 
know the Hubble length when it re-entered the horizon. 
It is useful to introduced a critical frequency, which corre- 
sponds to the mode re-entering the horizon when H£ — 1 . 
As discussed in Sec. Ill Bl this mode has e* = e c = \pl — 1. 
We can measure all other frequencies as a multiples of the 
critical frequency via 



I 

fc 



e*(V2- l){2 + e*)< 



1/4 



(A12) 



Here, we have used that e c a^ = e*a^. This is a cubic 
equation in e* that can be analytically inverted to give 
e* = e*(/// c )i which then yields H^£ as a function of 



f/fc- However, the general formula is complicated and 
not particularly enlightening. More interesting are the 
limits: 



( V2 + i)V3(£) 4/3 7 / > /c 



(A13) 



This yields that 
O 



= ;u ^4S-S 2 (m 2 (f) 



54.9, / < f c . 

H 2 [36.4(/// c ) 4 /3, f>f c . 

(A14) 

Hence, in order to understand the the frequency depen- 
dence of Slew, we need to know 1Z — IZ(f) from numeric 
simulations. 

Finally, we need to specify the actual value of f c . We 
make use of the fact that the universe expands adiabati- 
cally during radiation and matter domination. Therefore, 
conservation of entropy yields 



9s{T c )aiT* = gs(T )a 3 X 



(A15) 



Here, T c and To indicate the temperature at the critical 
epoch H c £ — 1 and today respectively, and gs measures 
the effective number of degrees of freedom in the mat- 
ter sector as a function of temperature. We can relate 
the temperature of radiation at the critical epoch to its 
density via 



Pc = Ae c 



6e c 



9cn 2 T ? 
30 



(A16) 



where we have written g c = gs{T c ). Then, it follows that 
fc H c a c 



fc 



i 1/3^ 

1 9o T o 



2ira 2vrao 2n£ g \^ Tc 

2(180^£ C )V4 J/12 (f) T °- (A17) 

To get a numerical answer, we can take Tq — 2.75 K, 
.go = 3.91 HI, and e c = -\/2- I- Then, 



f c = 3.3 x 10" 



n 1 \ 1/2 /mn\ 1/12 

0.1 mm \ ' / 100 x 



Hz. (A18) 

9c 
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